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Optimal Control of Transient Flow in Natural Gas Networks 

Anatoly Zlotnik' and Michael Chertkov^- and Scott Backhaus^ 


Abstract — We outline a new control system model for the dis¬ 
tributed dynamics of compressible gas flow through large-scale 
pipeline networks with time-varying injections, withdrawals, 
and control actions of compressors and regulators. The gas 
dynamics PDE equations over the pipelines, together with 
boundary conditions at junctions, are reduced using lumped 
elements to a sparse nonlinear ODE system expressed in vector- 
matrix form using graph theoretic notation. This system, which 
we call the reduced network flow (RNF) model, is a consistent 
discretization of the PDE equations for gas flow. The RNF 
forms the dynamic constraints for optimal control problems 
for pipeline systems with known time-varying withdrawals and 
injections and gas pressure limits throughout the network. 
The objectives include economic transient compression (ETC) 
and minimum load shedding (MLS), which involve minimizing 
compression costs or, if that is infeasible, minimizing the 
unfulfilled deliveries, respectively. These continuous functional 
optimization problems are approximated using the Legendre- 
Gauss-Lobatto (LGL) pseudospectral collocation scheme to 
yield a family of nonlinear programs, whose solutions ap¬ 
proach the optima with finer discretization. Simulation and 
optimization of time-varying scenarios on an example natural 
gas transmission network demonstrate the gains in security and 
efficiency over methods that assume steady-state behavior. 

I. Introduction 

New emissions restrictions and the resulting push towards 
cleaner electric power sources, as well as increased supplies 
of natural gas in the United States, have compelled the 
installation of gas-fired electric power plants for the vast 
majority of new generating capacity over the past 15 years 
[1]. Such generators can quickly adjust their output, hence 
they are often dispatched to balance out the fluctuating 
and highly variable production of uncontrollable renewable 
energy sources such as wind and solar [2], [3], This growing 
and increasingly time-variable natural gas consumption cre¬ 
ates a significant impact on the pressure and flow throughout 
associated natural gas transmission networks. These condi¬ 
tions contrast with historically slower and smaller variations 
in withdrawals by local distribution companies (LDCs), 
which allowed gas pipeline operators to assume, for day- 
ahead planning purposes, that consumption remains constant 
throughout the day. Today’s new complexities cause inter¬ 
actions on previously distinct spatiotemporal scales, which 
present risks and disruptive challenges that invalidate many 
traditional approaches for design, risk assessment, and oper¬ 
ation of natural gas transmission systems [4], [5], [6]. 
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Historically, natural gas was predominantly withdrawn 
from transmission systems by LDCs and industrial con¬ 
sumers in a predictable way and with relatively little variation 
over a day. It was traded using day-ahead contracts for 
fixed deliveries with the assumption that injections and with¬ 
drawals will remain nearly constant. Accordingly, early stud¬ 
ies [7], [8], [9] focused on optimizing steady-state gas flows, 
for which the state equations are algebraic relations. Recent 
efforts have improved and scaled up optimization techniques 
for similar problems [10], [11], [12], [13]. However, the 
steady-state assumption no longer represents realistic oper¬ 
ating conditions because of the economic and technological 
trends discussed above [14]. Intermittent dispatch of gas-fired 
power plants creates rapidly changing gas withdrawals from 
transmission pipelines, and results in stresses that are increas¬ 
ingly difficult to contain within system design limits using 
current ad-hoc methods. The growing reliance on natural gas 
for electricity production has led to strong coupling between 
electric power and natural gas infrastructures, and has created 
a need for secure gas transmission to prevent electric load 
shedding because of interrupted gas deliveries to generators 
[3]. This new challenge compels our investigation of new 
techniques for modeling and optimal control of dynamic 
flows of compressible gas in large-scale pipeline networks. 

The transient flow of natural gas in a transmission pipeline 
can be represented by simplifications to the Euler equations 
for compressible gas flow in one-dimension [15], [16]. For 
physically relevant pipeline parameters, this PDE system is 
defined over very large scales in both distance and time, and 
is highly nonlinear even with numerous physical modeling 
simplifications [17]. Transient flows in pipelines on the scale 
of thousands of miles are thus problematic to simulate, and 
many methods have been proposed [18], [19], [20]. Indeed, 
pipeline simulation is an area of active research of interest 
to gas system operators [21], 

The difficulty of characterizing gas network dynamics 
presents challenges for engineering, design, and operation of 
natural gas transmission systems under transient conditions. 
In the related optimization problems, the PDEs representing 
the dynamics are incorporated as constraints that must be 
satisfied over widely distributed space and time domains 
[22], and their nonlinearity makes computational tractabil- 
ity a challenge. Previous studies examined optimization of 
multi-day operations of gas pipeline networks involving 
transient flows [23], [24], [25], including work on economic 
model predictive control [26], [27]. In these studies the PDE 
constraints for gas flow are represented using implicit first- 
order schemes in space and time, which result in very large- 
scale problems because of the fine discretization required to 
adequately resolve transients on the time-scales of interest. 


In this manuscript, we develop a modeling and control 
framework that closely represents the physical phenomena 
in gas pipeline networks. We also present methods for sim¬ 
ulation and optimization of dynamic compressible gas flows 
over such systems using nodal actuators, which provides 
unprecedented gains in efficiency and scalability over the 
previous work listed above. For example, the optimal gas 
flow (OGF) [12] determines optimal compressor station set- 
points that balance constant injections and withdrawals over a 
network while satisfying system pressure constraints. Here, 
we extend this concept to the transient case. We examine 
the objectives of economic transient compression (ETC) and 
minimum load shedding (MLS), which involve minimizing 
compression costs or unfulfilled deliveries to customers with 
non-firm contracts, respectively. In the transient regime, these 
problems are formulated as PDE-constrained optimal control 
problems (OCPs). The PDE constraints are approximated 
by a new control system model, the reduced network flow 
(RNF), derived from a model reduction of gas network 
dynamics [19]. The approximation used for the RNF has 
been validated through comparison with an operator split- 
step numerical solution to the one-dimensional PDE system 
for a single pipe with transient boundary conditions [28]. 
The modified OCP is then approximated with a nonlinear 
program (NLP) using pseudospectral discretization [29]. The 
decision variables are coefficients of a polynomial expansion 
that approximates the solution to the OCP. Our approach 
provides several advantages relative to previous methods. 
The representation of continuous dynamics using polynomi¬ 
als gives spectral accuracy, yielding comparable fidelity using 
coarser discretization with far fewer decision variables. The 
RNF equations can also be integrated with an ODE solver 
to validate solutions to the NLP. 

The manuscript is organized as follows. In Section HU we 
summarize the physics of compressible gas flow in pipelines. 
Section [III] contains an optimal control formulation for dy¬ 
namic flows on networks subject to time-varying injections, 
withdrawals, and actions of compressors and regulators, and 
defines the ETC and MLS OCPs. In Section llVl we derive the 
RNF control system model, and provide consistency results. 
In Section [V] we summarize the Legendre-Gauss-Lobatto 
(LGL) pseudospectral collocation method for optimal con¬ 
trol. Section[Vl]describes implementation of the LGL scheme 
to approximate the reduced OCPs in Section [IV] with NLPs. 
Section I VIII contains a case study of the ETC and MLS 
problems for an example gas network. Section [YIIII contains 
a brief discussion of the results and their implications. 

II. Gas Pipeline Dynamics 


The dissipative flow of a compressible gas in a horizontal 
pipeline with slow transients that do not excite waves or 
shocks is adequately described by a simplification of the 
Euler equations in one dimension [18], [20], [16], given by 


dtp + d x p = 0 , 
d t p + a 2 d x p = - 


A p\p\ 


(1) 

( 2 ) 


The variables are mass flux p and density p, defined on a 
domain x £ [0, L] at time t. The parameters are the friction 
factor A, pipe diameter D, and speed of sound a. The term on 
the right hand side of ([2} aggregates friction effects. We have 
assumed that gas pressure p and density p satisfy the relation 
p = a 2 p with a 2 = ZRT, where Z, R, and T, are the 
gas compressibility factor, ideal gas constant, and constant 
temperature, respectively. Equation (0 is valid in the regime 
when changes in the boundary conditions are sufficiently 
slow to not excite propagation of sound waves. The gas 
dynamics on a pipeline segment are represented using 0- 
0 , with a unique solution when any two of the boundary 
conditions p(t, 0) = p 0 (t), p(t, 0) = p 0 (t), p(t,L) = p L {t ), 
or p(t,L) = pi,(t) are specified. 

Because of friction, the pressure of gas flowing through a 
pipeline gradually decreases, and is boosted by compressors 
so that it exceeds the minimum for delivery to customers. 
Compressor stations are controllable actuators used to ma¬ 
nipulate the state of the transmission system. The physical 
size of a station is very small compared to a pipeline, 
hence we model compressor/regulator action as conservation 
of flow and a multiplicative change in density at a point 
x = c. Specifically, p(t,c + ) = a(t)p(t,c~) and p(t,c + ) = 
p(t , c~), where a(t) is a time-dependent compression factor. 
We denote h(c ~) = lim x ^ c h(x) and h(c + ) = lim x x c h(x). 
The compression power is proportional to 

C oc p~ 1 \p(t, c)|(max{a(t), 1 } 2m — 1) (3) 

with 0 < to < (7 — l)/7 < 1 where 7 is the heat capacity 
ratio and p is the compressor efficiency [7], [12]. 

III. Network Flow Control Formulation 

We consider a network of pipelines that are connected 
at junctions where gas flow can be compressed, or gas can 
be withdrawn from or injected into the system. Between 
junctions, the mass flux and density evolve according to 
0 - 0 . This collection of segments connected at junctions 
is considered as a directed graph Q = (V, £), where each 
segment is an edge {i,j} £ £ in the set of edges £ that 
connects junctions i,j £ V in the set of nodes V. The 
instantaneous state within the edge {i,j} is characterized 
by the density pij and flux p, :l defined on a time interval 
T = [0, T\ and on the distance variable a Uj £ [0, L,;j ] = Cij, 
where L l3 is the length of edge Defining the domain 

of the PDE solution for the { 1 . j } 7 £ as 2?,, = Tx C t] , 
the state functions are expressed as ptj : T>ij —> K + and 
Pij : T>j —> R, where R + = (0, 00 ). We use a directed graph 
in order to denote for each edge a positive flow direction, 
i.e., if {i,j} £ £ then {j,i} £, which leads to the 

identity Pij(t,Xij) = —pji(t, Lij — ). We define the set 

of controllers C C £ x {+, —where {i,j} = {i,j, +} £ C 
denotes a controller located at node i £ V that augments the 
density of gas flowing into edge {i,j} £ £ in the + direction, 
while {j, i} = {i,j, —} £ C denotes a controller located at 
node j £ V that augments density into edge {i,j} £ £ in 
the negative direction. Compression is then modeled as a 
multiplicative ratio onj : T —> R+ for {i,j} £ C. 


2D p 





We denote by Sj : T —» R the density of gas entering the 
network from a node j £ Vs, where the set Vs denotes large 
supply terminals we call “slack” junctions, able to supply any 
mass flux at the given density. A mass flux withdrawal (or 
injection, if negative) at a junction j £ Vd = V \ Vs is 
denoted by dj : T —» R, where Vd is the set of demand 
(non-“slack”) nodes. The functions {&ij} {i,j}ec, Wj}jev D , 
and {sj jjev.s create nodal balance conditions of the form 


0) Djfc (j-')Pij (f, Lj ), 

Vj £ V D and {i,j}, {j , k} £ £, (4) 

dj(t)= Y Vijit’Lij)- Y Vv=M), VjGVn, (5) 
iev D kev D 

pij it, 0) = Si(t), V* £ Vs ( 6 ) 


To more conveniently represent the dynamics <[T]>- ( 0 for 
each edge, we apply the dimensional transformations 





V = 


V 

apo ’ 


(7) 


where l and po are nominal length and density, to yield the 
non-dimensional edge dynamics 


dt Pij T" dx^Pij — 
JJtPij V O x pij = 


0 , 

| '-p'l 'j | 

2Dij pij 


( 8 ) 

\/{i,j}££ (9) 


in which the hats have been omitted for readability. We 
henceforth use the non-dimensional units, and consider ®- 
<0 to represent flow dynamics on a pipeline segment. Ob¬ 
serve that this non-dimensionalization is not edge-dependent, 
hence the same factors £ and po are used system-wide. 

With this notation, we formulate two PDE-constrained 
OCPs for gas pipeline networks, for which the edge dy¬ 
namics ©-© and nodal conditions ©-© form dynamic 
constraints. Each is subject to transient withdrawals dj(t) 
for j £ Vd, available supply densities Sj(t) for j £ Vs, and 
box constraints on the density and compression of the form 


< Pij(t,Xij) < p” ax , 

V{i, j} £ £, 

(10) 

1 < ctij(t) < a^ ax , 

V {*> j} £ C- 

(11) 

For simplicity, we choose terminal conditions on the 
and control variables to be time-periodic, 

state 

Pij (0, Xij ) — Pij (T, Xij ), 

V{i, j} £ £, 

(12) 

(pij (0; Xij ) — (pij (V, Xij ) , 

V{i, j} £ £, 

(13) 

a io (0) = a.ij (T), 

V{*,j }£C, 

(14) 

though if the initial conditions are specified we may use 

__ fLij _ pLij 

Y / Pij{0,x)dx= Y 

{ij}e£ Jo VJ} ef 

/ Pij{T,x)dx, 

Jo 

(15) 


which is a periodic condition on total mass in the system. 

The ETC objective function, which aggregates compres¬ 
sion costs of the form © throughout the network, is 

J e =Y [ ((ma x{aij(t), l}) 2m - l)d£, (16) 

Vij 


where rpj is the efficiency of compressor {i , j} £ C. We 
formulate the ETC OCP for day-ahead operational planning, 
in which the system state is expected to repeat on a 24-hour 
period T. The decision functions are nodal controls 
for {i, j} £ C, and the problem is 

min J E in (fl 6 ] > 

s.t. dynamic constraints: ©-® 

nodal conditions: ©-© (17) 

density & control constraints: © -CB 
terminal constraints: m-m 

In the case that the ETC OCP ( fTTb does not have a feasible 
solution, we wish to determine the control protocol that 
comes closest to fulfilling the desired gas deliveries dj(t). In 
current practice, natural gas is purchased in day-ahead con¬ 
tracts for firm or non-firm delivery. Customers of the former 
type are guaranteed deliveries, while the latter may be cut off 
if their demand profile is thought to impact pipeline security, 
i.e. the density constraints (ITOt . We therefore formulate the 
MLS objective function in order to minimize the gap between 
the desired and achievable deliveries to non-firm customers 
at the nodes Ai = {ji, ■ ■ ■, jsj C Vd- Suppose that for 
j £ M, the desired withdrawal from the network is d*(t), 
and the MLS objective is given by 

J M = Y l c iWidjit) - d*(t)) 2 dt. (18) 

jeM 

Here, Jm is the sum of L 2 norms of delivery gaps weighted 
by the marginal quadratic costs Cj(t) of load-shedding at 
nodes j £ A4. Minimizing Jm leads to the desired distribu¬ 
tion of load-shedding across the pipeline system. 

The decision functions for the MLS OCP are nodal con¬ 
trols aij(t) for {*, j} £ C, as well as positive withdrawals 
dj(t) for j £ M, for which we must provide constraints 

dj( 0) = dj(T) and dj(t) > 0, Vj £ M. (19) 

The MLS problem is given by 

min Jm in (TTSl ) 

s.t. dynamic constraints: 

nodal conditions: © - © 

density & control constraints: © - © K ’ 
terminal constraints: ©-© 
delivery constraints: © 

The OCPs © and © are in general analytically in¬ 
tractable, and pose challenges for computational solution 
even for \£\ = 1 , i.e., a pipeline with no controllers or 
withdrawals except at the boundaries [30]. We proceed to 
derive the RNF control system based on a model reduction 
technique for flows on gas networks [19]. The RNF replaces 
the dynamic constraints in the OCPs (ITTb and (l20t . 

IV. Control System Model Reduction Of Gas 
Network Dynamics 

We construct a control system using a tractable yet accu¬ 
rate model of the network flow dynamics ©-© with nodal 
conditions ©-©, extending previous modeling work [19]. A 





Fig. 1. Relation E3 between nodal densities p^ and endpoint densities p9. 
and pfj, illustrated for a single edge (left) and for a joint (right). 


pipeline segment of (non-dimensional) length L is modeled 
as a lumped element by integrating (E])-© with respect to x, 
L 

(d t p + d x (p)dx = 0, (21) 

[ (d t p + d x p)dx = -2^- / ^Ux, (22) 

Jo 2D Jo P 

and then evaluating these integrals of <9 t , d x , and nonlinear 
terms using the trapezoid rule, the fundamental theorem of 
calculus, and averaging variables, respectively. This yields 


L 

2 

L 


(p° + p L ) = p°-p L , 
(p° + p L ) = p°-p L - 


\iL (gP + <p L )\ip° + ip L 


(23) 

(24) 


4D p° + p L 

where p°, gp and p L , tp L denote density and flux at the start 
and end of an edge, respectively, and p = ^p. 

to a network by defining input and 


ifiij and densities > 0 , p-j > 0 for 


We extend 
output flows 

each edge {ip} £ £. Then equations (I23l >- (l24b and nodal 
conditions ©-© then reduce to a differential-algebraic 
equation (DAE) system where the edge dynamics are 


Pit 


P% 




fij 


(25) 


<Pt 


■$3 


P% 


A ijt 
4A,- 


(vi 


+Vii)\ Vi 

Pii+Pi . 


+v°i\ 


(26) 


and the nodal constraints become 

n_ L 0 / Vj £ Vd, 

0 - atjkPij a jiPjk , | ^ ^ € (27) 

dj = Vij ~ H Vij> Vj £ Vd, ( 28) 

iGVd A:GVr> 

P% = Si, V i £ V S . (29) 

Here (127b represents continuity of pressure at junctions 
with jumps in the case of compression or regulation, (l28b 
represents flow balance at junctions, and (I25l i-(l26b represent 
flow dynamics on each segment. 

We apply the graph notation in Section|III]to express ( 125 l i- 
(l29b in matrix-vector form. Suppose that V = |V|, and assign 
to each node an index in [V], where [ N] = {1,..., N} for a 
positive integer N £ N. Define mappings 7 r° : £ —> [V] and 
n E : £ —» [V], which map the first and last vertices of an 
edge to the vertex ordering. Also suppose that E = \£\ and 
assign to each edge an index in [ 72 ], and define n e : £ —> \E\, 
which maps the edges to this ordering. We now assign to each 


node in the range of ~ v a unique internal density, and write 
a total nodal density state vector p N = [p^ ,..., Py) T ■ The 
mappings tt]! and n E can be used to state the dependence of 
p N on the variables { p ( k } and {pfj} and control variables, 

Pij ®ij TV) (ij) and Pi] P^T (if) ‘ (30) 

Equation (l30b will be used to state ( I25 | i-( |26| i in terms of 
nodal densities p N . We also write state vectors gp = 
(Vi,---,Vd) T and <P L = where gP k and 

p k are indexed by k = n e (ij). 

We then define the time-dependent weighted incidence 
matrix B : K s —>• R y by 

{ ctij edge k = n e (ij) enters node i, 

—aij edge k = n e (ij) leaves node i, (31) 
0 else 

as well as the incidence matrix A = sign (B). We define 
the collection of withdrawal fluxes d = (di,..., <7 m) T with 
M = |Vd|, where d k is negative if an injection. Also 
define the slack node densities as s = (si,..., sj,) T = 
{pf}j&Vs’ where b = |Vs|, and demand node densities as 
P = (pi, ■ ■ ■ ,Pm) t = (pf }jev D , so that b + M = V. 
Then let A S ,B S £ R &x£ denote the submatrices of rows 
of A and B corresponding to Vs, and let Ad,Bd £ R MxE 
correspond similarly to Vd- Next, let Al and Aq denote the 
positive and negative parts of Ad, so that Ad = Al + Aq. 
Also, define the diagonal matrices A, AT £ R ExE by A k k = 
L k and I\ k k = fA k/Dk, where L k , X k , and D k are the 
nondimensional length, friction coefficient, and diameter of 
edge k = n e (ij). Finally, define a function g : R e x Rf 
WL e by gj(x,y) = Xj\xj\/yj. Then (l25li-(l29b can be written 

d = A L <p L + A 0 (p o, (32) 

\B E \s+\B E \p=-Ak- l \{ip L -y 0 ), (33) 

\(J>l +<Po) = -A _1 (Sjs + Hjp) 

-Kg(X(g> L +ip 0 ),\Bj\s+\B E \p) (34) 

Furthermore, note that Ad = A k + Aq and | Ad\ = Al — Aq, 
so by defining <p = 2(<ps + <po) and <p- = \(ip L - <Po) 
we may replace ( l32b with d = Ad<p + \Ad\gr-, the right 
hand side of ([33} with —4A _1 (/?_, and the left hand side 
of ([34b with ip. Then multiplying (l33b by |Ad|A results in 
\A d \A\Bj\s + \A d \A\B E \p = -4\A d \p- = 4 {A d p-d). The 
DAE system d32t - (l34b may then be written as an ODE 

p = (I^IAIBJI)" 1 ^^ - d) - \Ad\A\Bj\s], (35) 

p = -A ~\Bj s + B E p) — Kg(ip, \B e \s + \Bj\p), (36) 

where p are nodal densities and p approximate mass flux on 
the edges. For a connected graph. Ad £ R MxE and Bd £ 
R MxE are full rank, and therefore |Ad|A|Hj| is invertible. 
Time-varying parameters are gas withdrawals d £ R M , input 
densities s £ R|j_, and compressions/regulations ctij £ C. We 
now give two consistency results. The RNF for a pipeline 
reduces to as the maximum spatial discretization step 

approaches zero, and reduces to the Weymouth equations 
[ 11 ] in the steady-state. 















Proposition 1: The RNF (|3D-(|3D is a consistent spatial 
discretization of the PDE ©-I® for a pipeline of dimen¬ 
sional length L, modeled as a chain of m segments of 
uniform length L/m, and which has no compressors or 
intermediate withdrawals (dj = 0 for j ^ m). 

Proof: The pipeline is represented as a graph Q with m = 
E, V = E + 1 , and each node at the points Xi = iL/m 
for i = 1 ,... ,m — 1 is connected to two edges and for 
* = 0 , 77 i is connected to one edge. It can be shown that 
W = |A d |A|f?J| is tridiagonal with W iti _i = L/(mt), 
Wi t i = 2L/(mt), and W^i+i = L/(mi), and matrix A d 
satisfies = 1, (A d ) iti+1 = -1. With pi(t) = p(t,Xi ) 

and = <p{t,Xi — \f) at time t, (f35l) - (l36l> yield 


f(Pi -1 + + Pi+l) — 


Al <Pi\<Pi\ 

2 D \{pi + Pi- 1 ) 


(37) 

(38) 


with scalar boundary conditions p± = s and <p m = d m - 
Taking m —> oo yields with the transformations ©. 


Remark 1: Proposition!]] implies that the solution to OD - 
( |3D using an implicit ODE integrator, which adapts the time- 
step to maintain stability and accuracy, will converge to the 
solution to I®-® as the space discretization M is increased 
[31]. Indeed, a simulation of the RNF approximation (|3P - 
(|3P for a standard pipeline model [20] with slow transients 
in input pressure and output flux yielded similar results as 
a solution of ©-® using an operator split-step method for 
hyperbolic PDE systems [28]. 

Proposition 2: When p = 0, f = 0, s = 0, d = 0, and 
aij = 0 for all {*, .)} €E C, equations (l35l >- (l36t reduce to the 
steady-state balance laws for a gas network [ 12 ]. 

Proof: Equation 1351) is reduced to A d p = d, which is 
nodal conservation of flow on demand nodes. Recall that 
p N contains all nodal densities, so equation (l36l > leads to 


AK<p © \p\ = (B T p N ) © (\B T \p N ), (39) 


where 0 denotes the point-wise vector product. Returning 
to {Pij} and {pfj} from the nodal densities using (l30t and 
reverting to the dimensional variables leads to the Weymouth 
equations for static gas networks [ 11 ], 

= <40) 

Remark 2: The potential equation < l40b . where {pT} for 
any slack node i £ Vs is given, together with the relation 
A d p = d, characterize the steady-state balance laws [12]. 

With the above results, we see that the RNF (l35l >- (l36l t can be 
used to represent both the PDE model ([l}-® for a dynamic 
network and the static equilibrium equations. 

To implement the RNF in the transient ETC (fTTI ) and MLS 
(l20l > OCPs, we rewrite the constraints and objective functions 
in Section[III]in terms of the nodal densities p and edge fluxes 
ip used in (l35l >-(l36l>. The inequality constraint (ITOt becomes 

pf n < MOM*) < P“ ax , (41) 


and the time-periodic terminal conditions (fT2l> - (fT3l> become 

p(0)=p(T), p(0) = <p(T), (42) 

and total mass conservation (fl5T> becomes 
l T A(|77j|(s(0) - s(T)) + \Bj\(p(0) - p(T))) = 0 (43) 


where 1 is a vector with one in each coordinate. The ETC 
objective function (fl~ 6 l> becomes 

Je='^2, [ ^ ((maxla^t), l}) 2m -l)df (44) 

while the objective (fT 8 l> for MLS remains unchanged. Using 
the RNF, the reduced ETC OCP takes the form 

min Je in (l44l > 

s.t. RNF constraints: (l35l > — (l36t 

density & control constraints: (ED. CD 
terminal constraints: <E3.([H 


The reduced MLS OCP is given as 

min Jm in CD 

s.t. RNF constraints: (l35l ) — (l36l> 

density & control constraints: (ED. CD (46) 
terminal constraints: ED-d 
delivery constraints: CD 

Although we have reduced the instantaneous states to 
the vectors p(t) and ip(t) from the continuous functions 
Pij{t,Xij) and <pij(t,Xij) used in Section!]]]] the OCPs (ED 
and (ED require optimization on the space of functions ( t ) 
for all {i,j} £ C. Therefore, we employ a method for time- 
discretization of such problems to NLPs. 


V. PSEUDOSPECTRAL OPTIMAL CONTROL 

We review a numerical scheme for transcribing OCPs into 
NLPs using pseudospectral discretization. Consider the OCP 

min J(x,u) = I C(t,x(t),u(t))dt, (47) 

u ' Jo 

s.t. x(t) = f(t,x(t),u(t)), (48) 

e(®(0),x(r)) = 0, (49) 

g(x(t),u(t)) < 0, (50) 

on T = [0, T]. Here C £ C K is in the space C K of 
continuous functions with k classical derivatives, and the 
dynamic constraints / £ C^~ x are in the space ~ 1 of 
n-vector valued C K_1 functions, with respect to the state, 
x(t) £ K", and control, u(t) £ R m . The functions e and 
g are terminal and path constraints, respectively. Finally, the 
admissible set for controls u includes the piecewise Cf func¬ 
tions on T. For details, refer to [29]. We now derive a direct 
collocation procedure for constructing a finite-dimensional 
nonlinear program that approximates the problem (I47l >- (|5()1 >. 

We use Lagrange polynomial interpolation to approximate 
the states x and controls u at a set of collocation points tk- 

x{t) « x N (t) = Y^k=o (51) 

u{t) w ujv(f) = Lfclo (52) 





Lagrange interpolation polynomials satisfy lk(U) = Ski, 
where Ski is the Kronecker delta function [32]. It follows 
that x{t k ) = x N (t k ) = x k and u{t k ) = u N (t k ) = u k , 
so the physical meaning of the interpolating polynomial 
coefficients x k and u k are the values of the state and 
control variables at the collocation points. Those points are 
chosen so that the integral in (l47l > and the derivative in 
(l48l > are computed accurately. The former is approximated 
using Legendre-Gauss quadrature, leading to the choice of 
Legendre polynomials as the orthogonal basis. Furthermore, 
the Legendre-Gauss-Lobatto (LGL) quadrature points are 
chosen to include endpoints of the interval, permitting the 
terminal constraints to be specified within the scheme. The 
LGL quadrature rule for a function / : [— 1,1] —> R. is 

1 N pi 

/ f(t)dt Wi= ii{t)dt, (53) 

J-i i=1 J -1 

and is exact if / G P 2 Ar~i and the nodes t t G y lgl , where 
P 2 Ar-i denotes the set of polynomials of degree at most 
2TV — 1 and where Y LGL = {t , : L^{ti) = 0, i = 1,... TV — 
1} lj{—1,1} are the TV + 1 LGL nodes determined by the 
derivative of the .Y ,h order Legendre polynomial, L^{t), and 
the interval endpoints [32], Because the OCP is defined on 
T = [0, T], whereas Legendre polynomials form a basis on 
[—1,1], we re-scale time by t = (2 1 — T)/T. 


We re-write the Lagrange interpolating polynomials on the 
LGL collocation nodes in terms of the Legendre polynomial 
basis, to provide the scheme with derivative and spectral 
accuracy properties of orthogonal polynomials. Given tk G 
Y lgl , we can express the Lagrange polynomials as [33] 


f (f \ = i 

U N(N + l)L N (t k ) t-t k 


(54) 


The derivative of (BTt at t 3 G Y LGL is then 


d_ 

dt 


N N 

^ ^ XkXk (tj ) — ^ ( DjkXk, 
fc= 0 fc= 0 


(55) 


where D is the constant differentiation matrix with elements 
D ik = 4(fi). Using (ED, d52), (ED, and 03, the OCP 03- 
(ES) is transcribed as the following NLP, in which the de¬ 
cision variables are the polynomial interpolation coefficient 
vectors x = (xq, • ■ •, xjv ) and u = (uq, ..., un)- 


N T 

min J(x,u) = ^ —C(x k ,u k )w k 


(56) 

k—0 



N T 

S.t. ^ ^ DikXk — " 9 "^ 'U'i) ^ ^ — 0? 1? ■ ■ 

..,7V 

(57) 

k—0 



e(x 0 ,x N ) = 0 


(58) 

g(x k ,u k ) <0 V k = 0,1,., 

..,7V 

(59) 


The solutions to (I56ll -(l59l> converge to extrema of (I47ll - (l5()l) 
as TV —> oo at an exponential rate because of the spectral 
accuracy of polynomial approximations [29]. 


VI. Implementation 

The model reduction in Section [IV] is in essence a spatial 
discretization of the the pipeline flow dynamics in Section 
m on a graph Q = (V, £), which incorporates boundary 
conditions at network nodes. For the model to accurately ap¬ 
proximate the true PDE dynamics, this discretization must be 
sufficiently fine. Thus to translate from the PDE-constrained 
problems ( fl7l) or (l20l ) to the reduced OCPs (l45l > and (l46l t. 
we first create a modified graph Q = (V,£) by adding 
nodes such that all edges of £. are shorter than a maximum 
length l, which can also be used as the non-dimensional 
constant. It has been observed [34], and we have confirmed 
[28], that l = 10 km is sufficient to adequately represent 
transients of interest for typical transmission pipelines. The 
graph Q is then used to create the RNF (I35li-(l36l>. which is 
used for © or (l46l >. The system state variables are then 
p G R M and ip G where M = \Vd\ and E = \£\, so 
that x = ( p , ip) is the state vector of interest. The control 
variables are for {i,j} G C, where C = \C\ is the 
number of compressors/regulators. In the case of the MLS 
problem, the withdrawals d :l for j G M are variables as 
well, and £ = \M.\ are the number of non-firm loads where 
shedding may occur. The control vector u contains all a, :l 
and dj which are to be determined by the system operator. 
The OCPs (1451) and (l46l) can then be expressed in the form 
of (07])-®, and then approximated as an NLP using the 
procedure described in Section[V] In the NLP (156| )-( |59| ). each 
element of the vector-valued functions x and u is expressed 
using N +1 Lagrange interpolation coefficients, which leads 
to (M + E + C ) x (IV+ 1) and (M + E + C + Z) x (IV+ 1) 
variables for the ETC and MLS OCPs, respectively. 

To guarantee a smooth, physically relevant solution, we 
add a penalty on the square of the L 2 norms of derivatives 
of the compression ratios to the objective function: 

( T 

Js{a) = p ll«vll 2 = M / (“iiW) 2 ^ 

{i,j}e£ Jo 

N / N 

E E Dmk^ijh 

{i,j}e£m=0 \k —0 

where ctij = {ctijo, • • •, oeijN) T are interpolation coefficients 
for compression function {i,j} G C, and // is a relative 
weight, for which an effective empirical value is N. The 
cost term Jg is eliminated when discussing objective values. 

The optimal control scheme is implemented computa¬ 
tionally as follows. First, all system parameters including 
network structure and constraints as well as interpolation 
coefficients of time-varying withdrawals and injections are 
used to build MATLAB functions for the objective, con¬ 
straints, and their gradients with respect to the decision 
variables. These are provided, along with random initial 
conditions that satisfy inequality constraints, to the interior- 
point solver IPOPT version 3.11.8 running with the sparse 
linear solver ma57 [35]. Convergence of optimization for 
the case studies below requires only minutes because the 
gradients are provided to the solver, and the constraint 






Fig. 2. Example network (not to scale). Numbering indicates nodes (blue, 
above/right), edges (black, below/left), and compressors (red, below/right). 
Thick and thinner lines indicate 36” and 25” pipes, respectively. Pressure is 
bounded between 500 and 800 psi on all pipes. Friction factor and sound 
speed are A = 0.01 and a = 377.968 m/s. 

Jacobian provided to ma57 has under 3% non-zero entries. 
VII. Examples 

Current industry practice is to assign compression set- 
points in an ad-hoc manner with the stipulation that gas with¬ 
drawals are constant throughout the day [36]. Because actual 
gas flows may be up to 80% above or below the planned- 
for rates, set-points must be chosen very conservatively, so 
that pressures may drop far below the rated minimum, which 
leads to load-shedding and gas price spikes. 

Solving the ETC problem addresses such issues by ac¬ 
counting for transient withdrawals, which are usually known 
by gas system operators on a day-ahead basis [37]. By 
utilizing time-dependent dynamical information, the transient 
ETC is much more likely to have a solution that is actually 
valid for operations, and load-shedding will be unnecessary. 
An intermediate formulation, which we call “quasi-static”, 
is examined for comparison, where constant compression 
set-points are chosen to satisfy pressure constraints given 
dynamic withdrawals. When the transient ETC problem does 
not have a feasible solution, the transient MLS problem can 
determine the most efficient protocol for load-shedding. The 
latter leads to the ultimate utilization of pipeline network 
capacity while ensuring that pressure limits are not exceeded. 

We consider an example system in Figure [2] consisting of 
a tree with 25 nodes connected by 24 edges with a total 
length of 477 km, containing C = 5 compressors, and a 
single slack node j = 1. For more accurate representation of 
transients, artificial nodes are added so that pipe segments 
have maximum length 10km, resulting in V = 62 nodes 
(yielding M = 61 non-slack nodes), and E = 61 edges. Gas 
is supplied at the minimum pressure of 500 psi at the slack 
node, to be immediately compressed into the network. 

ETC Optimization. The compression ratios for the system 
in Figure [2] are optimized to solve the ETC problem ( 145| i for 
the withdrawal profiles in Figure [3f. Two slowly-changing 
profiles on node sets {12,19} and {6, 8} represent residential 
and industrial use, while profiles for the sets {18,24} and 
{25} represent common single-cycle gas turbine operations. 
Node 13 has a constant injection. Thus the slow and fast 
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Fig. 3. ETC Optimization, (a) Transient compression solutions (color), with 
scaled extremal system pressure (thick dashed line) and limits (dotted line); 
(b) Quasi-static solution - maximum pressure relaxed by 25%; (c) Static 
solution - pressure drop as seen in practice; (d) Pressures at nodes (MPa); 
(e) Fluxes on edges (kg/m 2 /s); (f) Withdrawal/injection profiles (kg/m 2 /s). 

demand profiles account for 34.5% and 65.5% of the to¬ 
tal consumption, respectively. The problem is solved with 
N = 25 time points, so the total number of variables is 
(M + E + C) x (N + 1) = 3302, and solution takes 
under 5 minutes on a laptop computer. Figure [3 }l shows the 
transient compression solution, as well as resulting scaled 
maximum and minimum pressures, which fall within the 
bounds as indicated. Figures [3}> and [3J; show these results 
for the quasi-static solution and the fully static solution, 
respectively. A feasible quasi-static solution is obtained only 
when the maximum pressure constraint is relaxed with a 
25% increase. The scaled extremal pressures that result 
when these solutions are applied to simulations with the 
transient withdrawals dramatically violate the desired limits, 
as observed in current pipeline operations [36]. 

MLS Optimization. Consider the same scenario as for the 
ETC problem, except the loads at non-firm customers at 
nodes {18,24} are increased so that there is no feasible 
solution. Therefore the MLS objective is used, which adds 
(N + 1) x E = 52 optimization variables. A priority 
weighting of ci§(f) = 024 (f) = 1 is used. The desired and 
maximal non-firm deliveries are shown in Figure [4{>, and 
the associated compression solutions are given in Figure [4jt. 
With the right control protocol, the system utilization can 
be significantly increased even over the delivery profiles in 
Figure [3f and still yield a feasible solution. 
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Fig. 4. MLS Optimization, (a) Compression solutions (color), scaled ex¬ 
tremal system pressure (thick dashed line), and scaled pressure limits (dotted 
line); (b) desired (solid) and maximal non-firm (dashed) deliveries. 

VIII. Conclusion 

We have developed a new framework for modeling and 
optimal control of compressible gas flow through pipeline 
networks with time-varying injections, withdrawals, and con¬ 
trol actions of compressors and regulators. The inclusion of 
information about transient parameters into a physically rep¬ 
resentative model and an efficient and tractable optimization 
scheme are shown to facilitate a dramatic and unprecedented 
improvement in both the capacity and security of pipeline 
operations with respect to current industry practice. The 
economic transient compression (ETC) formulation provides 
a cost efficient solution when the desired mass transfer is 
feasible, and the minimum load shedding (MLS) objective 
minimizes unfulfilled deliveries if they cannot all be met. 
Moreover, the solutions produced by the optimization scheme 
are validated by direct simulation of the control system 
model. In addition, our technique leverages the inherent 
sparsity of the problem for efficient scaling to larger systems. 
Implementation in practice would dramatically increase the 
effective capacity of gas pipeline systems, and save signifi¬ 
cant resources now used to provide unutilized margins. 
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